home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byte1286.arc / AITKEN.F77 < prev    next >
Encoding:
Text File  |  1986-12-05  |  10.0 KB  |  304 lines
  1.       PROGRAM GAUSS
  2. C
  3. C  FORTRAN-77 PROGRAM TO APPROXIMATE DEFINITE INTEGRALS.  VERSION: 4-14-86
  4. C
  5. C  DAVID M. SMITH
  6. C  MATHEMATICS DEPARTMENT
  7. C  LOYOLA MARYMOUNT UNIVERSITY
  8. C  LOS ANGELES, CA  90045
  9. C
  10. C
  11.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  12.       DIMENSION TABLE(20)
  13. C
  14. C             UNIT NUMBERS:  KW - SCREEN OUTPUT
  15. C                            KR - KEYBOARD INPUT
  16. C                            KF - OUTPUT TO FILE 'GAUSSQ.OUT'
  17. C
  18.       KW = 6
  19.       KR = 5
  20.       KF = 8
  21.       OPEN(KF,FILE='GAUSSQ.OUT',STATUS='NEW')
  22. C
  23. C             FORTRAN-77 COMPILERS ALLOW THESE 5 INPUT VALUES TO BE
  24. C             ENTERED IN FREE FORMAT.  I.E., TO ENTER NFCT = 2 TYPE THE
  25. C             2 IN COLUMN 1 AND THEN TYPE CARRIAGE RETURN.  TO ENTER
  26. C             A = 1 TYPE 1. OR 1.0 THEN RETURN.  THE DECIMAL POINT
  27. C             SHOULD BE TYPED FOR THE REAL VALUES (A AND B).
  28. C
  29. C             MOST FORTRAN-66 COMPILERS WHICH SUPPORT IF-THEN-ELSE
  30. C             REQUIRE THAT THE INTEGERS BE TYPED RIGHT-JUSTIFIED,
  31. C             SO THAT TO ENTER NFCT = 2 TYPE 4 BLANKS AND THEN THE 2
  32. C             IN COLUMN 5.
  33. C
  34.   110 WRITE (KW,120)
  35.   120 FORMAT(/' Enter (1 item per line):  NFCT / KL / IOFLAG / A / B.'
  36.      *       /7X,'(Format for NFCT,KL,IOFLAG is I5, for A,B is F25.0.'/
  37.      *       /' Integrate function NFCT from A to B  (NFCT=0 to quit)',
  38.      *       /' There will be KL lines in the Gauss',
  39.      *       ' column of the table.'/
  40.      *       ' IOFLAG = 1 causes the output to be saved in file',
  41.      *       ' GAUSSQ.OUT.'/)
  42. C
  43.       READ (KR,130) NFCT
  44.   130 FORMAT(I5)
  45.       IF (NFCT.EQ.0) STOP
  46.       READ (KR,140) KL,IOFLAG,A,B
  47.   140 FORMAT(I5/I5/F25.0/F25.0)
  48. C
  49.       IF (NFCT.LE.0) STOP
  50. C
  51.       KPRT = 1
  52.       IF (IOFLAG.EQ.1) KPRT = 2
  53.       CALL GAUSSQ(A,B,NFCT,KL,TABLE,KPRT,KF,KW)
  54.       DO 170 J = 1, 20
  55.          WRITE (KW,150)
  56.   150    FORMAT(/' ENTER 1 FOR THE NEXT AITKEN COLUMN, 0 TO QUIT.',
  57.      *          '  (I1)')
  58.          READ (KR,160) KOPT
  59.   160    FORMAT(I1)
  60.          IF (KOPT.NE.1) GO TO 110
  61.          CALL AITKEN(TABLE,KL,KPRT,KF,KW)
  62.          KL = KL - 2
  63.          IF (KL.LT.3) GO TO 110
  64.   170    CONTINUE
  65.       GO TO 110
  66. C
  67.       END
  68.       SUBROUTINE GAUSSQ(A,B,NFCT,KL,TABLE,KPRT,KF,KW)
  69.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  70. C                                                DAVID M. SMITH  4-14-86
  71. C
  72. C  ITERATIVE GAUSS QUADRATURE.  FUNCTION NUMBER NFCT IS INTEGRATED FROM
  73. C  A TO B PRODUCING KL LINES IN TABLE USING
  74. C  1, 2, 4, 8, 16, ..., 2**(KL-1) SUBINTERVALS.
  75. C
  76. C  IF KPRT IS 1 THE TABLE IS WRITTEN ON UNIT KW.
  77. C  IF KPRT IS 2 THE TABLE IS ALSO WRITTEN ON UNIT KF.
  78. C
  79. C  IN THE TABLE WHICH IS WRITTEN:
  80. C
  81. C  N        IS THE NUMBER OF SUBINTERVALS USED,
  82. C  GAUSS    IS THE GAUSS QUADRATURE APPROXIMATION TO THE INTEGRAL,
  83. C  TOTAL NF IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS DONE SO FAR,
  84. C  EST S.D. IS A CONSERVATIVE ESTIMATE OF THE NUMBER OF SIGNIFICANT
  85. C           DIGITS WHICH ARE CORRECT IN THAT LINE
  86. C
  87. C  THE EXTERNAL FUNCTION CALLED IS F(X,NFCT) FOR THE FUNCTION BEING
  88. C  INTEGRATED.
  89. C
  90.       DIMENSION TABLE(20)
  91. C
  92. C             SET EPS TO THE SIZE OF THE ROUNDING ERRORS FOR A GIVEN
  93. C             MACHINE.  EPS = 1.0E-7 FOR 32 BIT MACHINES (SINGLE
  94. C             PRECISION) OR 1.0E-16 (DOUBLE PRECISION).
  95. C
  96.       EPS = 1.0E-16
  97. C
  98. C             CHECK FOR ERRORS.
  99. C
  100.       IF (KL.LT.1) THEN
  101.           WRITE (KW,110) KL
  102.   110     FORMAT(/' ERROR IN GAUSSQ.  KL=',I5,'.  IT SHOULD BE AT ',
  103.      *           'LEAST 1.')
  104.           DO 120 J = 1, 20
  105.   120        TABLE(J) = -3.1E+31
  106.           RETURN
  107.       ENDIF
  108. C
  109. C             DO GAUSS QUADRATURE FOR THE INTEGRAL OF F(X) FROM A TO B.
  110. C
  111.       NLINES = KL
  112.       SQRT3 = DSQRT(3.0D0)
  113.       NSUBS = 1
  114.       NFUNCT = 0
  115. C
  116.       IF (KPRT.GE.1) THEN
  117.           WRITE (KW,130) NFCT,NLINES
  118.           IF (KPRT.GE.2) WRITE (KF,130) NFCT,NLINES
  119.   130     FORMAT(/' GAUSS INTEGRATION OF FUNCTION ',I3,'.  THERE',
  120.      *           ' ARE',I3,' LINES IN THE TABLE.')
  121.           WRITE (KW,140) A,B
  122.           IF (KPRT.GE.2) WRITE (KF,140) A,B
  123.   140     FORMAT(/' THE LIMITS OF INTEGRATION ARE:'/3X,D23.16,
  124.      *           '  TO  ',D23.16/)
  125.       ENDIF
  126. C
  127.       DO 180 JLINE = 1, NLINES
  128. C
  129. C             THE KTH SUBINTERVAL IS (AK,BK).  THE GAUSS APPROXIMATION
  130. C             FOR THE INTEGRAL IS:  XH*(F(XM-XR) + F(XM+XR)), WHERE
  131. C             XH = (BK-AK)/2,   XM = (AK+BK)/2,   XR = XH/(2*SQRT(3)),
  132. C             AK = A + (K-1)*(B-A)/NSUBS,   BK = A + K*(B-A)/NSUBS.
  133. C
  134.          XH = (B-A)/NSUBS
  135.          XH2 = XH/2.0
  136.          XR = XH2/SQRT3
  137.          START1 = A - XH2 - XR
  138.          START2 = A - XH2 + XR
  139.          SUM = 0.0
  140. C
  141.          DO 150 K = 1, NSUBS
  142. C
  143.   150       SUM = SUM + F(START1+K*XH,NFCT) + F(START2+K*XH,NFCT)
  144. C
  145.          SUM = SUM*XH2
  146. C
  147.          NFUNCT = NFUNCT + 2*NSUBS
  148.          IF (KPRT.GE.1) THEN
  149.              IF (JLINE.LE.1) THEN
  150.                  WRITE (KW,160) NSUBS,SUM,NFUNCT
  151.                  IF (KPRT.GE.2) WRITE (KF,160) NSUBS,SUM,NFUNCT
  152.   160            FORMAT(' N =',I6,'  GAUSS=',D23.16,'     TOTAL NF=',I5)
  153.              ELSE
  154.                  RELERR = EPS
  155.                  IF (SUM.NE.0.0) RELERR = DABS(TABLE(JLINE-1)-SUM)
  156.      *                                    / DABS(SUM)
  157.                  IF (SUM.EQ.0.0 .AND. TABLE(JLINE-1).NE.0.0) RELERR =
  158.      *               DABS(TABLE(JLINE-1)-SUM) / DABS(TABLE(JLINE-1))
  159.                  IF (RELERR.LE.0.0) RELERR = EPS
  160.                  NUMSD = -DLOG10(RELERR)
  161.                  IF (NUMSD.LT.0) NUMSD = 0
  162.                  WRITE (KW,170) NSUBS,SUM,NFUNCT,NUMSD
  163.                  IF (KPRT.GE.2) WRITE (KF,170) NSUBS,SUM,NFUNCT,NUMSD
  164.   170            FORMAT(' N =',I6,'  GAUSS=',D23.16,'     TOTAL NF=',I5,
  165.      *                  '     EST S.D.=',I3)
  166.              ENDIF
  167.          ENDIF
  168. C
  169.          NSUBS = 2*NSUBS
  170.          TABLE(JLINE) = SUM
  171.   180    CONTINUE
  172. C
  173.       RETURN
  174.       END
  175.       SUBROUTINE AITKEN(TABLE,KL,KPRT,KF,KW)
  176.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  177. C                                                DAVID M. SMITH  4-14-86
  178. C
  179. C  AITKEN EXTRAPOLATION.  THE VALUES TABLE(1) TO TABLE(KL) ARE USED FOR
  180. C  THE EXTRAPOLATION AND THE RESULTS ARE RETURNED IN TABLE(1) TO
  181. C  TABLE(KL-2).  KL MUST BE AT LEAST 3.
  182. C
  183. C  IF KPRT IS 1 THEN THE KL-2 VALUES ARE WRITTEN ON UNIT KW.
  184. C  IF KPRT IS 2 THE TABLE IS ALSO WRITTEN ON UNIT KF.
  185. C
  186. C  IN THE TABLE WHICH IS WRITTEN:
  187. C
  188. C  LINE     IS THE LINE NUMBER,
  189. C  AITKEN   IS THE AITKEN EXTRAPOLATION APPROXIMATION TO THE INTEGRAL,
  190. C           USING 3 VALUES FROM THE PREVIOUS COLUMN.  FOR EXAMPLE,
  191. C           LINE 1 OF AN AITKEN COLUMN USES LINES 1-3 OF THE PREVIOUS
  192. C           COLUMN, LINE 4 USES LINES 4-6, ETC.
  193. C  EST S.D. IS A CONSERVATIVE ESTIMATE OF THE NUMBER OF SIGNIFICANT
  194. C           DIGITS WHICH ARE CORRECT IN THAT LINE
  195. C
  196.       DIMENSION TABLE(20)
  197. C
  198. C             SET EPS TO THE SIZE OF THE ROUNDING ERRORS FOR A GIVEN
  199. C             MACHINE.  EPS = 1.0E-7 FOR 32 BIT MACHINES (SINGLE
  200. C             PRECISION) OR 1.0E-16 (DOUBLE PRECISION).
  201. C
  202.       EPS = 1.0E-16
  203. C
  204.       IF (KL.LT.3) THEN
  205.           WRITE (KW,110) KL
  206.   110     FORMAT(//' ERROR IN AITKEN.  KL=',I5,' IS LESS THAN 3.'/)
  207.           RETURN
  208.       ENDIF
  209. C
  210.       IF (KPRT.GE.1) THEN
  211.           KLM2 = KL - 2
  212.           WRITE (KW,120) KLM2
  213.           IF (KPRT.GE.2) THEN
  214.               IF (KLM2.GT.1) THEN
  215.                   WRITE (KF,120) KLM2
  216.   120             FORMAT(/' AITKEN COLUMN.',I4,' ENTRIES.'/)
  217.               ELSE
  218.                   WRITE (KF,130) KLM2
  219.   130             FORMAT(/' AITKEN COLUMN.',I4,' ENTRY.'/)
  220.               ENDIF
  221.           ENDIF
  222.       ENDIF
  223. C
  224.       KLM1 = KL - 1
  225.       DO 160 JLINE = 2, KLM1
  226.          TOP = (TABLE(JLINE+1) - TABLE(JLINE))**2
  227.          BOT = TABLE(JLINE+1) - 2.0*TABLE(JLINE) + TABLE(JLINE-1)
  228. C
  229.          IF (BOT.EQ.0.0) THEN
  230.              TABLE(JLINE-1) = TABLE(JLINE+1)
  231.          ELSE
  232.              TABLE(JLINE-1) = TABLE(JLINE+1) - TOP/BOT
  233.          ENDIF
  234. C
  235.          IF (KPRT.GE.1) THEN
  236.              JLM1 = JLINE - 1
  237.              IF (JLINE.LE.2) THEN
  238.                  WRITE (KW,140) JLM1,TABLE(JLM1)
  239.                  IF (KPRT.GE.2) WRITE (KF,140) JLM1,TABLE(JLM1)
  240.   140            FORMAT(' LINE',I3,'  AITKEN=',D23.16)
  241.              ELSE
  242.                  RELERR = EPS
  243.                  IF (TABLE(JLM1).NE.0.0) RELERR = DABS(TABLE(JLM1-1) -
  244.      *                                   TABLE(JLM1))/DABS(TABLE(JLM1))
  245.                  IF (TABLE(JLM1).EQ.0.0 .AND. TABLE(JLM1-1).NE.0.0)
  246.      *               RELERR = DABS(TABLE(JLM1-1)-TABLE(JLM1)) /
  247.      *                        DABS(TABLE(JLM1-1))
  248.                  IF (RELERR.LE.0.0) RELERR = EPS
  249.                  NUMSD = -DLOG10(RELERR)
  250.                  IF (NUMSD.LT.0) NUMSD = 0
  251.                  WRITE (KW,150) JLM1,TABLE(JLM1),NUMSD
  252.                  IF (KPRT.GE.2) WRITE (KF,150) JLM1,TABLE(JLM1),NUMSD
  253.   150            FORMAT(' LINE',I3,'  AITKEN=',D23.16,
  254.      *                  '     EST S.D.=',I3)
  255.              ENDIF
  256.          ENDIF
  257. C
  258.   160    CONTINUE
  259. C
  260.       RETURN
  261.       END
  262.       FUNCTION F(X,NFCT)
  263. C
  264. C  COMPUTE F(X) FOR FUNCTION NUMBER NFCT.
  265. C  THIS CALLING SEQUENCE ALLOWS MANY DIFFERENT FUNCTIONS TO BE USED
  266. C  EASILY DURING ONE RUN.
  267. C
  268.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  269. C
  270. C             IF NFCT IS OUTSIDE THE RANGE OF CURRENTLY DEFINED
  271. C             FUNCTIONS THEN RETURN THE FUNCTION VALUE -9.9E+20.
  272. C             ANY OUTPUT TABLES CONSISTING ENTIRELY OF THIS VALUE ARE
  273. C             ALMOST CERTAINLY CAUSED BY SENDING THE WRONG VALUE OF
  274. C             NFCT.
  275. C
  276.       F = -9.9E+20
  277.       IF (NFCT.LT.1 .OR. NFCT.GT.6) RETURN
  278. C
  279.       GO TO (1,2,3,4,5,6),NFCT
  280. C
  281.     1 F = DSQRT(DEXP(X)-1.0D0)/DSIN(X)
  282.       RETURN
  283. C
  284.     2 F = DABS(X-0.3)**0.33333
  285.       RETURN
  286. C
  287. C             THIS FUNCTION IS 1/(X**4 + X**2 + 1)
  288. C             THE FORM USED BELOW EXECUTES FASTER.
  289. C
  290.     3 F = 1.0/((X*X+1)*X*X + 1)
  291.       RETURN
  292. C
  293.     4 T = X*X + 1.0
  294.       F = T/(T*X*X + 1)
  295.       RETURN
  296. C
  297.     5 F = DSQRT(X)/(X - 1.0) - 1.0/DLOG(X)
  298.       RETURN
  299. C
  300.     6 F = 2.0*X*X/((X-1.0)*(X+1.0)) - X/DLOG(X)
  301.       RETURN
  302. C
  303.       END
  304.